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The properties of a front between two different phases in the presence of a smoothly inhomogeneous 
external field that takes its critical value at the crossing point is analyzed. Two generic scenarios are 
studied. In the first, the system admits a bistable solution and the external field governs the rate in 
which one phase invades the other. The second mechanism corresponds to a second order transition 
that, in the case of reactive systems, takes the form of a transcritical bifurcation at the crossing 
point. We solve for the front shape and its response to external white noise, showing that static 
properties and also some of the dynamics features cannot distinguish between the two scenarios. 
The only reliable indicator turns out to be the fluctuation statistics. These take a Gaussian form in 
the bifurcation case and a double-peak shape in a bistable system. The results of a recent analysis 
of the morphogenesis process in Drosophila embryos are reanalyzed and we show, in contrast with 
the interpretation suggested by Krotov, et al. [T], that the plausible underlying dynamics is bistable 
and not bifurcational. 

PACS numbers: 87.10.Mn,87.23.Cc,64.60.Ht,05.40.Ca 


The theory of first and second order phase transitions 
is a well-established part of statistical physics, and its 
generalization to out-of-equilibrium problems, like glassy 
behavior and percolation-like transitions, has also re¬ 
ceived a lot of attention during the last decades. How¬ 
ever, most of these works focused on the case where 
the control variable (like temperature) is homogenous 
in space. Only recently has the equilibrium properties 
of thermodynamic systems with a spatially varying tem¬ 
perature that takes its critical value only in a localized 
region begun to be studied [2]. 

Coincidentally, an out-of-equilibrium process with the 
same spatial characteristics was suggested as the under¬ 
lying mechanism behind one of the most fundamental 
and universal aspects of developmental biology, embry¬ 
onic morphogenesis. Measuring the expression levels of 
gap genes along the anterior-posterior axis of Drosophila 
embryos, Krotov, et al. ^ revealed a strong negative cor¬ 
relation between spatial expression levels of gene pairs, 
with an increase in the expression level of gene A along 
the spatial axis associated with a decrease in expression 
of gene B. Krotov, et al. interpret their findings as ev¬ 
idence for competition between A and B; i.e., they as¬ 
sume that the transcription activity of A is reduced when 
gene B is active and vice versa. In the simplest two-gene 
model, the primary maternal morphogens enhance the 
production of A to the left of a crossing point and the 
production of B to its right. At the crossing point, where 
there is no preference to any of the genes, the system is 
at “criticality”, showing larger fluctuations with a fast 
manifold associated with conservation of the global ac¬ 
tivity and a slow manifold that corresponds to the differ¬ 
ences between concentrations of A and B. To model this 
process, Krotov, et al. implemented a simple reaction- 
diffusion model with a spatial gradient of the control pa¬ 
rameter (maternal morphogens) and white noise, arguing 
that the analytical and numerical results obtained from 


this model fit the empirically observed data. 

However, we believe there is a loophole in the logic of 
Krotov, et al. Their model, which involved two coupled 
logistic populations with a transcritical bifurcation at the 
crossing point is indeed one generic possibility, but there 
is an alternative generic scenario, wherein the system is 
bistable and the external held reaches the stall point at 
the transition. In this letter we analyze these two sce¬ 
narios and solve for the front width and the fluctuation 
spectrum at the crossing point. Doing that, one realizes 
that the second scenario, the bistable front model, bet¬ 
ter captures the features of the morphogenetic process 
considered in [1]. 

Qualitatively speaking, the bistable scenario is the ana¬ 
log of a simple hrst order transition system, although (as 
we shall see below) one should make a distinction between 
its equilibrium and out-of-equilibrium properties. Imag¬ 
ine a water-ice mixture in three dimensions, say, where 
the temperature depends on x and T{x = 0) = T^, where 
Tm is the melting temperature. In the right half-space 
water invades ice and in the left ice grains grow in wa¬ 
ter, so the ice-water front will be trapped around x = 0. 
Due to thermal or other fluctuations the front will move 
back and forth in a region determined by Xt, leading to 
a characteristic fluctuation spectrum at a: = 0; below 
we will distinguish between the instantaneous shape of 
the front and its average over time and will analyze the 
distribution of fluctuations. 

The bifurcation model, which is the one considered in 
[T], is slightly more complicated. This model is a nonequi¬ 
librium continuous (“second order”) transition, with a 
transcritical bifurcation at a; = 0. To imitate the gene 
expression case one needs two competing fields. A generic 
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set of PDEs that yields the required behavior is 
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Where a and b are the expression levels of the corre¬ 
sponding genes A and B, C{x) is the background field 
(maternal morphogenes) that switches sign at zero and C 
is, say, white noise. 

In the absence of noise, C = 0, in the regime C{x) > 0 
the fixed point (FP) of this system is a = 0, 6=1 
while if C{x) < 0 the FP is a = 1, 6 = 0. The diffusion 
term, however, couples the left and the right regions, and 
the expression level must be a smooth function of x that 
approaches the FPs at large \x\ and takes the symmetric 
value o = 6= l/2at the crossing point. 

An appropriate and quite generic choice of the external 
field profile is C{x) = tanh(a;/a;t), where Xt sets the scale 
of the external field gradient. Plugging 6 = 1 — a into 
Eqs. Q one finds in the deterministic limit (from here 
on we normalize K to unity), a = D'V^a — C{x)a{l — a); 
expanding C aX x Xt the equation for the stable front 
is then 


a''{y)-ya{y)[l-a{y)]=0, ( 2 ) 

where y = x{Dxt)~^^^ and spatial derivatives are taken 
with respect to y. 

Eq. (U) describes a deterministic, spatial voter model 
with selection [3], where the value of y determines the 
preference towards one of the ’’alleles” (opinions). The 
front width is proportional to [DxtY/^] shape and 
the fluctuation spectrum at the front will be determined 
below. 

It appears to be instructive to contrast the bifurcation 
model ([^ with the other generic scenario, a bistable sys¬ 
tem. Let us consider the simplest model of a bistable 
system with a crossing point, described by a spatially 
inhomogeneous Ginzburg-Landau (GL) Equation: 

^(x, t) = D\7^cj){x, t) + 4>{x, t)[l — (/)(x, t)\[(j){x, t) — C{x)]. 

( 3 ) 

When C is a constant, cf) admits three spatially homoge¬ 
nous FPs, two stable FPs at (^ = 0 and (p = 1 and an 
unstable FP at G, the (j) = 1 invades the zero phase if 
C < 1/2 and zero invades if C > 1/2. C = 1/2 is 
the melting point, or the stall point of the GL front; 
at the melting point, when the system evolves from in¬ 
homogeneous initial conditions like (/ = 0 for x < 0 and 
(/ = 1 for X > 0, it relaxes to the stable front solution, 
(j)o = [tanh(x/-\/8D) -I- l]/2. Accordingly, even when C 
is X dependent, as in the case C{x) = tanh(x/xt) con¬ 
sidered above, as long as the intrinsic width of the front 
\/%D is much smaller then Xt, the shape and the width of 
the front will be essentially independent of the external 



FIG. 1: The static properties of the deterministic fronts: 
panel (a) shows both static front profiles, ao(x) (pink, wider) 
and 4’o{x) (black line), for D = 0.4 and xt = 0.5. Clearly 
there is no essential difference between the two. In panel (b) 
the width of the bifurcation front was plotted against 
(for fixed Xt = 1, results are open circles, red line is a linear 
fit) and in panel (d) the same quantity is plotted against xt 
for fixed D — 10“^; as predicted, the bistable front width is 
linear in the square root of D and in independent of xt for 
when Xt is larger than the natural width of the front. Panels 
(d) and(e) show the same relations for a bifurcation front {D 
varies for xt = 1, Xt varies for D = 10“^), demonstrating the 
{DxtY^^ scaling. 


field (see Figure [^. Although both scenarios support a 
stable front, the dependence of its width on the exter¬ 
nal parameters is different: in a bifurcation system this 
width scales like and depends on the width of the 
crossing region xt, while in a bistable system the width 
scales like and is independent of Xt when xt is large. 

However, when a front is observed, as in the experi¬ 
ment discussed by [T], and one would like to determine 
the underlying mechanism, the utility of diagnostic tools 
based on static properties of the front is quite limited. 
In experiments it is quite difficult to change D oi Xt - to 
do that, the dynamics on the molecular level should be 
manipulated - so one cannot measure the dependency of 
the front width on D. Worse than that, it turns out that 
the front shape is almost identical in both cases. In the 
bifurcation model = 1/2-1- x/Vl6D close to x = 0, 
meaning that the front satisfies, to first order in x and </>, 
UV^(/(x) -I- {x/y/\&D)(j){x)[l — (/(x)] = 0, i.e., it is equiv¬ 
alent, up to a constant, to the bifurcation front solution 
([^, so the differences between </o(a:) and the solution of 
Eq. (§ (denoted hereon as ao(x)) are extremely small, 
as demonstrated in Figure]^ Without measuring the dif¬ 
fusion constant of the underlying morphogenesis molecu¬ 
lar agents, or monitoring the front profile to a very high 
degree of accuracy, one cannot use static properties to 
distinguish between the two possible scenarios. 

Not only the static properties of the front are practi¬ 
cally useless as an indicator of the underlying mechanism. 
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the same is true for some dynamical aspects of the dy¬ 
namics. In [I], for example, the location of the crossing 
point was identified (assuming an underlying bifurcation 
scenario) by a peak in the anticorrelations between the 
densities of a and b, and the existence of a slow and fast 
manifold was demonstrated by a scatter plot of the fluc¬ 
tuations, showing that the sum a -I- & is kept almost fixed 
through time but the differences a — b fluctuate strongly. 
Indeed, the same features are also a characteristic of a 
bistable scenario. To show that, we have developed a 
two-species model that supports bistability (in the one 
species case ^ the features demonstrated in [T] are em¬ 
barrassingly trivial, since the field (j) should be interpreted 
such that a = (j) and b = 1 — (j), so the sum is fixed and 
the anticorrelations are guaranteed in advance). 

To construct a simple two species bistable model we 
define S{x,t) = a{x,t) + b{x,t) and Q{x,t) = a{x,t) — 
b{x,t), and the local dynamics satisfies 

S = S{a-S) Q = {Q-C{x))iS^-Q^) (4) 

so the stable FPs correspond to S = a and Q = ±a and 
for constant C there is an unstable FP at Q = C. Both 
reactants a and b diffuse with a diffusion constant D. As 
shown before, the bistable front has an intrinsic width 
and its shape is independent of the external field param¬ 
eter xt as long as Xt ^ VSD. Accordingly we simulate 
this system with antiperiodic boundary conditions and 
without an external field. The correlation function of a 
and b, together with a scatter plot of the fluctuations at 
the crossing point, are depicted in Figurej^ showing that 
both models have very similar qualitative behavior. 

However, the sharp-eyed observer may notice a subtle 
qualitative difference between the scatter plots of fluctu¬ 
ation amplitudes. In the bifurcation model simulations 
the points appear to have higher density in the middle 
(around [0.5,0.5], which is the steady state value of the 
front at the crossing point), while the simulation of the 
bifurcation model yields higher concentration of fluctua¬ 
tion points close to the two extremes a = 0, 6 = 1 and 
0 = 1, 6 = 0. This is not an accident, and provides a 
crucial hint: the two mechanisms, bifurcation and bista¬ 
bility, admit qualitatively different fluctuation statistics. 
In the bistable scenario, due to the absence of the exter¬ 
nal gradient, the noise causes the front to move back and 
forth freely around the crossing point, so at a; = 0 the sys¬ 
tem is almost always either it the (j) = 1 state or the (j) = 0 
state, leading to a fluctuation spectrum with two peaks 
at zero and one and a dip at 1/2. The bifurcation mecha¬ 
nism, on the other hand, yields only a single peak around 
the steady state value ao{x = 0) = b{x = 0) = 1/2. 

To quantify this, we consider first the fluctuations 
around the steady state front of the bifurcation model, 
ao(x) + S(x, t) in the presence of an external white noise. 
Linearizing Eq. © to the first order in 6, and taking 
into account the front shape close to the crossing point, 
ao(a;) ~ l/2-|-cia;/(Dxt)^/^, where ci is an 0{1) constant. 



FIG. 2: The correlation between two timeseries, a(x, t) and 
h{x, t), for a fixed distance x from the crossing point, is plotted 
against x for both bifurcation model (full blue line) and the 
bistable model (red line with full circles). In both models the 
anticorrelations between fluctuations of a and b are peaked 
at x = 0 (main panel). In the inset, a scatter plot of the 
instantaneous fields at different times is presented in the a — 
b plane (without fluctuations there should be one point at 
1/2,1/2). Points from the bistable model are represented by 
filled red circles, bifurcation model points are empty and blue. 
Fast (a-|-6) and slow (a—b) manifolds are clearly seen. Results 
were obtained with xt = 10, D = 1/2 and noise amplitude 
0.04. 

one obtains a dynamical equation for the fluctuations of 
the bifurcation model: 

S{x, t) = D\7^6{x, t) — Kx'^S{x, t) + ({x, t) (5) 

where k = and C is a white noise, 

C(x,t) = 0 and ({x,t)(^{x',t') = Ad{x — x')6(t — f) where 
an overbar represent an average over realizations. Ex¬ 
panding S in terms of normalized quantum harmonic os¬ 
cillator wavefunctions, S{x,t) = Pm{t)'<Pra{x), and us¬ 
ing their orthonormality properties one obtains, 

/3m(t) = + rjii) i (6) 

with Tm = 2'/DK{m + 1/2) and r]{t) is, again, white 
noise. Every coefficient pm is subject to an Ornstein- 
Uhlenbeck process and its probability distribution func¬ 
tion is given by a Gaussian with zero mean and variance 
A/Fm- An immediate result is that S itself is a zero mean 
Gaussian, i.e., that the fluctuation density histogram is a 
Gaussian centered at ao(x = 0) = 1/2. Indeed one can do 
even better and calculate the variance of this Gaussian, 

Var{S) = 4’mi0)Var{rm) (7) 

m even 

_ A /H y- {(m-iyAp _ Ar(5/4) 

V D m!(m+l/ 2 ) r(3/4)cf« V D ' 

In a bistable system the situation is completely differ¬ 
ent. As explained above, as long as Xt is significantly 
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FIG. 3: Fluctuations statistics; a histogram (unnormalized) 
of a{t) values at the crossing point for the bifurcation (red 
line with filled circles) and bistable (blue) models. In both 
cases noise leads to deviations from the steady state value 
a = 0.5, however in the bifurcation case these deviations are 
distributed normally around the average while the bistable 
system distribution is bimodal. Simulation parameters are 
identical to those specified in the caption of Figure 


larger than the internal width of the front, one can re¬ 
place the external field (with exponentially small correc¬ 
tions in a finite system) by antiparallel boundary con¬ 
ditions at ±oo, and the the fluctuations admit a zero 
(Goldstone) mode since the location of the front is trans- 
lationaly invariant. Accordingly, one finds the crossing 
point either in the a phase or in the b phase, with fluctu¬ 
ations due to the effect of noise on any of these phases. 
As a result the histogram of fluctuations amplitude, in¬ 
stead of being a Gaussian around 1/2, has two peaks that 
correspond to the two attractive fixed points of the ho¬ 
mogenous problem. These features are demonstrated in 
Figure where the strong qualitative difference, allow¬ 
ing for easy discrimination between the two scenarios, is 
manifest. On the other hand, when xt is much smaller 
than the natural width of the front, even the bistable sys¬ 
tem loses its translational invariance property, the front 
is trapped by the external field and cannot change signif¬ 
icantly its spatial location, and the resulting fluctuation 
spectrum is peaked at 1/2. In such a case we cannot offer 
a simple method to distinguish between the two alterna¬ 
tives mechanisms. 

Turning back to the work of Krotov, et al. [T], the 
results from the experiment they analyzed clearly show 
a crossing regime with anticorrelations between a and b 
with fast and slow manifold, however, as we explained 
here, this cannot reveal the nature of the dynamics gov¬ 
erns the system. The only simple qualitative indicator is 
the histogram of the amplitude of the fluctuations, and 
their results (Fig 3G of [1]) clearly show a double peak 
structure, meaning that the underlying dynamics is evi¬ 
dently bistable, equivalent to a first order transition with 
an external field (primary maternal morphogenes) that 
changes the “temperature” such that the melting tem¬ 


perature marks the crossing point. This appears to rule 
out the bifurcational interpretation suggested in [1]. 

Finally, in the context of the bistable model we would 
like to stress the difference between two possible defini¬ 
tions of a front separating two phases. The analysis fol¬ 
lowed Eq. ([^ regards the instantaneous shape of a front, 
i.e., the typical shape of a snapshot of the crossover re¬ 
gion. In contrary one may define the time averaged, or 
the ’’equilibrium” front, wherein that the a density, say, 
is averaged at x over a long time span and the resulting 
front is the profile of {a{x)), where (...) represents the 
time, or equivalently the thermodynamic, average. 

The width of an equilibrium front under smoothly 
varying external field was analyzed by [2] in the context 
of a 2d q-state Potts model, these authors found that for 
q > 4 when the system has a first-order transition, the 
width of the front {a{x)) scales as a;/ . To understand 
and generalize their result, let us consider an equilibrium 
system at the transition point. Starting from a homoge¬ 
nous A state, B phase droplets with the same bulk energy 
are nucleated and shrink only due to the surface tension. 
As a result, the larger the droplet, the more stable it is; 
monitoring the phase at a certain point x for long time 
one finds {a{x)) = 1/2, independent of the location of the 
measurement. Considering a randomly moving front, like 
the one described above, one arrives at the same conclu¬ 
sion. 

What limits the size of these droplets, hence deter¬ 
mining the width of the equilibrium front, is the external 
field gradient. If phase A invades the region x > 0 (where 
phase B is prefered) by a compact semisphere of radius 
A, the bulk energy cost U of such a droplet is, 

pR d-l ^ ^ 

U(x — x^) ^ —dx ^ -, (8) 

Jq 

meaning that the width of the equilibrium front scales 
like (the scaling when d = 2, was found in 

0)- Accordingly, the width of the equilibrium front does 
depend on xt, as opposed to the instantaneous front. In 
any case, the hallmark of a bistable system is the double 
peak of the fluctuation spectrum, not the shape of the 
front. 

The problem considered here, a front pinned by smooth 
spatial gradient of an external field, appears to be quite 
generic. Beyond the experimental results that were con¬ 
sidered in ^ , it may be relevant to the effects of environ¬ 
mental gradient on the genetic heterozygosity of a pop¬ 
ulation (see, e.g., 0) and on the species richness, gene 
transfer and speciation in ecological communities along 
such a gradient (known as an ecotone or ecoline) [S1|B]. In 
particular, the distinction between a stable, bifurcational 
front and the wandering front characterizing a bistable 
scenario may be very relevant to the rate of gene flow 
and to the chance of ecotonal species to survive. Further 
studies of these phenomena, and in particular an appro¬ 
priate classification of front dynamics using fluctuation 
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statistics, may shed a new light on many fundamental 
processes both in physics and in the life sciences. 
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